%% ==================== [n-D estimate of the output statistic] =========================
clear
close all
polOrd=3;  
dimension=2;
inVarMat=n1nDpolyRoots(polOrd,dimension);
outVarMat=112+3*(inVarMat(:,1)*2).^2+exp(inVarMat(:,2))-0.2*inVarMat(:,1).^3 ...
+0.5*sin(inVarMat(:,2));
samplePointNum=1000; %sampling points
                                                                     
[PCmean,PCvar,MC_PCmean,MC_PCvar,MC_PCdata,MCinput]=GetStatistics ...
    (inVarMat,outVarMat,polOrd,samplePointNum);
syms x1 x2
 
MCoutput=112+3*(x1*2).^2+exp(x2)-0.2*x1^3+0.5*sin(x2);
for i=1:2
    MCoutput=subs(MCoutput,sym(['x',num2str(i)]),MCinput(:,i)); 
end
MCoutput=double(vpa(MCoutput,4));
[Xpdf,X] = ksdensity(MC_PCdata); 
MCmean=mean(MCoutput);
MCvar=var(MCoutput);
[X1pdf,X1] = ksdensity(MCoutput); % find the probability density function of the random variable
%% ========================== [plot] =================================
plot(X1,X1pdf,'LineWidth',1)
xlim([100,200])
hold on
plot(X,Xpdf,'---','LineWidth',1);
title('Probability density of random variables')
xlabel('Random variable Y')
ylabel('probability density \it{f}')
legend('MCpdf','PCpdf')